Daremos por hecho que podemos generar variables aleatorias independientes e identicamente distribuidas con distribución Uniforme$(0,1)$.
Empezamos por la distribución más sencilla posible: $\text{Bernoulli}(p=0.5)$.
Teniendo acceso a una variable aleatoria Uniforme$(0,1)$; ¿cómo simularían la variable aleatoria $X\sim \text{Bernoulli}(0.5)$?
En efecto $X=\mathbb{1}_{\left\{ U \leq 0.5 \right\}}$ se distribuye $\text{Bernoulli}(p=0.5)$
En general si tenemos una distribución discreta $X$ tal que para $x_i\in \left\{ x_1,x_2,\ldots ,x_n\right\}$
\begin{align*} \mathbb{P}\left[X=x_i \right]=p_i \end{align*}con
\begin{align*} \sum_{i=1}^n p_i =1. \end{align*}¿Cómo se puede simular $X$ haciendo uso de una variable aleatoria $U\sim\text{Uniforme}(0,1)$?
En efecto podemos definir
\begin{align*} X=\begin{cases} x_1 & \text{ si } U\leq p_1 \\ x_2 & \text{ si } p_1 < U \leq p_1 + p_2 \\ \vdots \\ x_n & \text{ si } p_1 + p_2 + \ldots + p_{n-1}< U \end{cases} \end{align*}La metodología anterior puede ser generalizada haciendo uso de $F^{\leftarrow}$, la inversa generalizada de la función cumulativa de distribución $F$ asociada a la variable alaeatoria $X$ que se quiere simular.
Recordamos que la función de distribución cumulativa de una variable aleatoria $X$ está dada por
\begin{align*} F(x) = \mathbb{P}\left[ X \leq x \right]. \end{align*}Una función de distribución cumulativa debe ser no-decreciente, continua por la derecha y tal que
\begin{align*} \lim_{x\to-\infty}F(x)=0 \end{align*}y \begin{align*} \lim_{x\to \infty}F(x)=1. \end{align*}
Definimos la inversa generalizada de una función cumulativa de distribución como
\begin{align*} F^{\leftarrow}(u)=\inf\left\{ x \, : \, F(x) \geq u \right\} =\min\left\{ x \, : \, F(x) \geq u \right\}, \quad 0<u<1. \end{align*}Dada $F$ una función de distribución cumulativa:
Si $X\sim \text{Exponencial}(\lambda) \implies F(x)=1-e^{-\lambda x}\implies F^{\leftarrow}(x)=F^{-1}(x)=-\frac{\log(1-x)}{\lambda}$, en pracica se usa $X=-\frac{\log(U)}{\lambda}$ para simular variables aleatorias exponenciales.
using Plots, RandomNumbers.Xorshifts
r = Xoroshiro128Plus(0x1234567890abcdef)
Xoroshiro128Plus(0xe7eb72d97b4beac6, 0x9b86d56534ba1f9e)
n = 1000
x = -log.(rand(r,n));
plot(p)
plot(p)
Si se condiciona $X$ a $X\in (a,b)$ para $a<b$, entonces para $x\leq b$
\begin{align*} \mathbb{P}[X < x | X\in(a,b)] = \frac{ \mathbb{P}[a < X < \min\left\{ x , b\right\} ] }{\mathbb{P}[a<X<b]} = \frac{ F(x) -F(a)}{F(b)-F(a)}. \end{align*}Se sigue que $F^{-1}(F(a)+U(F(b)-F(a)))\sim X|X\in (a,b)$.
Sea $X$ una variable aleatoria continua, con función de densidad de probabilidad $f(x)$, de la cual deseamos generar simulaciones. La idea del método de aceptación-rechaza es usar una variable aleatoria auxiliar $Y$, de la cual podemos simular y un criterio de aceptación de la variable aleatoria tal que la distribución de $Y$ condicionada al criterio de aceptación coincida con la distribución de interés asociada a $X$.
En particular consideraremos $Y$ una variable aleatoria auxiliar continua con densidad de probabilidad $g$ y una constante $0<M\in \mathbb{R}$ tales que $f(x) \leq Mg(x) \, \forall x\in \mathbb{R}$.
1) Simule $Y\sim \mathcal{L}(g)$ y $U\sim \text{Uniforme}(0,1)$.
2) Acepte $X=Y$ si se cumple que $$U\leq \frac{f(Y)}{Mg(Y)};$$ en caso contrario regrese a 1).
Observemos que una función de densidad de probabilidad $f$ siempre se puede escribir como
\begin{align*} f(x) = \int_0^{f(x)}\mathrm{d}u \end{align*}por lo que $f$ es la densidad marginal de la distribución conjunta
\begin{align*} (X.U) \sim \text{Uniforme}\left\{ (x,u) \, : \, 0<u<f(x) \right\} \end{align*}es equivalente a $$ (X,U)\sim \text{Uniforme}\left\{ (x,u) \, : \, 0<u<f(x) \right\} $$
Simular una dupla $(X,U)$ como en el teorema fundamental es sencillo si simulamos $X\sim \mathcal{L}(f)$ y $U | X \sim \text{Uniforme}(0,f(X))$, sin embargo esto no es de utilidad cuando la simulación de interes corresponde a $\mathcal{L})(f)$. Alternativamente, una metodología de aceptación-rechazo interesante en este contexto consiste en considerar $(Y,V)\sim \text{Uniforme}(M)$ con $M$ un conjunto tal que $\left\{ (x,u) \, : \, 0<u<f(x) \right\}\subseteq M$ y aceptar si la dupla $(Y,V)\in \left\{ (x,u) \, : \, 0<u<f(x) \right\}$,
Para ilustrar tal método de aceptación-rechazo consideramos una densidad $f$ tal que $$ \int_a^b f(x)\mathrm{d}x=1 $$ para $a<b$ números reales. Entonces podemos considerar $m=\max\left\{ f(x) \, : \, a \leq x \leq b \right\} $ y $M=\left\{ (y,v) \, : \, 0<v<m \wedge a<y<b \right\}$. Usamos el siguiente algoritmo de simulación.
plot(p)
plot(p)
Anteriormente consideramos el conjunto $M=\left\{ (y,v) \, : \, 0<v<m \wedge a<y<b \right\}$ que nos ayudo a definir variables auxiliares distribuidas uniformemente sobre $M$. Más en general podemos considerar.
$M=\left\{ (y,v) \, : \, 0<v<m(y) \right\}$. Donde sin pérdida de generalidad podemos escribir $m(y)=Mh(y)$ para $0<M<\infty$ y $h$ una densidad de probabilidad.
Observemos que con $V | Y \sim \text{Uniforme}(0,m(Y))$ la regla de aceptación "$ X=Y \iff V\leq f(Y)$" tiene sentido para un algoritmo de simulación de $f$ si $f(x)\leq m(x)\, \forall x\in \mathbb{R}$; aún más en tal caso se tiene que $\mathcal{L}(X)=\mathcal{L}(f) \iff h=g$.
De la observación anterior para $g$ densidad de probabilidad y $M$ constante positiva tales que $f(x)\leq Mg(x) \, \forall x\in \mathbb{R}$ obtendríamos el siguiente algoritmo
1) Simule $Y\sim \mathcal{L}(g)$ y $V\sim \text{Uniforme}(0,Mg(y))$.
2) Acepte $X=Y$ si se cumple que $$V\leq f(Y)$$ en caso contrario regrese a 1).
Que coincide con el algoritmo de aceptación-rechazo al considerar que $U\sim \text{Uniforme}(0,1)\iff KU\sim \text{Uniforme}(0,K).$
plot(p)
Para el algoritmo de aceptación-rechazo es necesario poder evaluar $Mg(x)$ y deseable que esta evaluación tenga un bajo costo computacional. Más aún el algoritmo puede ser modificado para basarse en la evaluación de $cf(x)$ para alguna constante $c>0$, esto permite hacer simulación cuando la constante de normalización de la densidad $f$ es inaccessible o costosa de evaluar.
1) Simule $Y\sim \mathcal{L}(g)$ y $U\sim \text{Uniforme}(0,1)$.
2) Acepte $X=Y$ si se cumple que $$U\leq \frac{cf(Y)}{Mg(Y)};$$ en caso contrario regrese a 1).
En este caso la probabilidad de aceptar es $\frac{c}{M}$ y se mantiene que $X\sim \mathcal{L}(f)$.
Podemos considerar
\begin{align*} f(x) & \propto e^{-x^2/2}\left( sin(6x)^2 +3cos(x)^2 sin(4x)^2 + 1 \right) \\ g(x) & \propto e^{-x^2/2} \end{align*}y $M=5$
plot(p)
p=histogram(y_vec[accept_ind], nbins=250, normalize=true, label="")
title!("Histograma de simulaciones de f con aceptación rechazo")
plot(p)
Si podemos evaluar de manera no costosa una función $g_l$ tal que
$$ g_l(x) \leq f(x) \leq M g_m(x) $$con $M>0$ constante positiva y $g_m$ densidad de probabilidad tales que se puede implementar el algoritmo de aceptación rechazo. Entonce tiene sentido proponer el siguiente algoritmo cuando evaluar $g_l$ es menos costoso computacionalmente que evaluar $f(x)$.
1) Simule $Y\sim \mathcal{L}(g_m)$ y $U\sim \text{Uniforme}(0,1)$.
2) Acepte $X=Y$ si se cumple que $$U\leq \frac{g_l(Y)}{Mg_m(Y)};$$ en caso contrario proceda a 3).
3) Acepte $X=Y$ si se cumple que $$U\leq \frac{f(Y)}{Mg_m(Y)};$$ en caso contrario regrese a 1).
Decimos que una función de de densidad de probabilidad es log-cóncava si $h(x)=log(f(x))$ es cóncava, si $f$ es suficientemente suave esto es equivalente a que $h'(x)=\frac{\mathrm{d}h}{\mathrm{d}x}$ sea estrictamente decreciente o $h''(x)=\frac{\mathrm{d}^2 h}{\mathrm{d}x^2}<0$. Para este tipo de densidades es sencillo construir funciones que "ensandwichen" a la densidad de interés.
plot(p)
plot(p)
plot(p)
Recordamos que si $h$ es diferenciable la recta tangente a la gráfica de $h$ en $x_i$ está dado por $T(x)=h(x_i)+h'(x_i)(x-x_i)$.
Dados $x_i<x_{i+1}$, las rectas tangentes a $h$ en $x_i$ y $x_{i+1}$ se intersectan en el punto dado por
\begin{align*} z_i = \frac{h(x_{i+1}) -h(x_i) -x_{i+1}h'(x_{i+1}) +x_ih'(x_i) }{h'(x_i)-h'(x_{i+1})} \end{align*}Por lo que dados puntos $T_n = \left\{ x_1, \ldots , x_n \right\}$ definimos para $x\in(z_{i-1},z_i]$
$$ M g_u(x) = h(x_i) + h'(x_i)(x-x_i) $$con $i\in \left\{1, \ldots, n\right\}$, $z_0$ el mínimo del soporte de $h$, en caso de que $h$ tenga soporte acotado por debajo ó $-\infty$ en caso contrario, y $z_k$ es el máximo del soporte de $h$ en caso de que $h$ tenga soporte acotado por arriba ó $\infty$ en caso contrario.
Por otro lado definimos para $x\in [x_i,x_{i+1}]$
$$ g_l(x) = \frac{(x_{j+1}-x)h(x_j)+(x-x_j)h(x_{j+1})}{x_{j+1}-x_j} $$y si $x<x_1$ o $x>x_n$ definimos $g_l(x)=-\infty$.
1) Incialize $n$ y $T_n$ (el conjunto puntos nudo).
2) Genere $X\sim \mathcal{L}\left( g_u(x) \right)$ y $U \sim \text{Uniforme}(0,1)$.
3) Si $$U\leq \frac{ g_l(X) }{Mg_u(X)}$$ acepte $X$; en caso contrario si $$U\leq \frac{ f(X) }{Mg_u(X)} $$ acepte $X$ y actualize $n=n+1$, $T_{n+1}=T_n \cup \left\{ X \right\}$; en caso contrario regrese a 2).
@time prueb_prev = adaptive_accept_reject_vprev([-1.1,1.1],10000,h,dh);
29.384700 seconds (1.61 G allocations: 26.664 GiB, 14.47% gc time)
@time prueb = adaptive_accept_reject([-1.1,1.1],10000,h,dh);
0.606691 seconds (12.94 M allocations: 319.221 MiB, 12.00% gc time)
@time randn(10000);
0.000094 seconds (6 allocations: 78.359 KiB)
plot(p)
plot(p)
plot(p)
La distribución gaussiana inversa generalizada está determinada por la siguiente función de densidad de probabilidad
\begin{align*} f_{GIG}(x)= \frac{\left(a/b\right)^{p/2}}{2\mathcal{K}_p\left(\sqrt{ab}\right)} x^{p-1}e^{-\left(ax +b/x \right)/2}\mathbb{1}_{\left\{ x >0 \right\}} \end{align*}con $a,b>0$, $p\in \mathbb{R}$ y $\mathcal{K}_p$ la función modificada de Bessel del segundo tipo.
Observemos que
\begin{align*} f_{GIG}(x) \propto x^{p-1}e^{-\left(ax +b/x \right)/2}\mathbb{1}_{\left\{ x >0 \right\}} \end{align*}por lo que existe $c>0$ tal que \begin{align*} h(x)=\log \left( c f_{GIG}(x) \right) =(p-1)\log(x)-\frac{ax}{2} - \frac{b}{2x} \end{align*} satisface \begin{align*} h''(x)=-\frac{(p-1)}{x^2} -\frac{b}{x^3} \end{align*} por lo que $h$ es cóncava para $p\geq1$ con $x$ restringida a $(0,\infty)$.
plot(p)
plot(p)
Regresión Poisson
Supongamos que tenemos muestras $(Y_1,X_1), \ldots , (Y_n,X_n)$ tales que
\begin{align*} Y_i | a,b,X_i & \sim \text{Poisson}\left( e^{aX_i + b} \right) \\ a & \sim \text{Normal}(0,\sigma^2) \\ b & \sim \text{Normal}(0,\tau^2) \end{align*}Si consideramos observaciones $Y_i=y_i$, la distribución posterior de $a,b|\pmb{Y}=\pmb{y},\pmb{X}=\pmb{x}$ entonces tiene densidad dada por
\begin{align*} \pi(a,b|\pmb{y},\pmb{x}) &\propto \pi(\pmb{y}|\pmb{x},a,b) \pi(a) \pi(b) \\ & \propto exp\left( a\sum_{i=1}^n y_i + b\sum_{i=1}^n y_i x_i - e^{a}\sum_{i=1}^n e^{x_i b} -\frac{a^2}{2\sigma^2} -\frac{b^2}{2\tau^2} \right) \end{align*}Entonces
\begin{align*} log\left( \pi(a|\pmb{y},\pmb{x},b) \right) & = a\sum_{i=1}^n y_i - e^{a}\sum_{i=1}^n e^{x_i b} -\frac{a^2}{2\sigma^2} \\ log\left( \pi(b|\pmb{y},\pmb{x},a) \right) & = b\sum_{i=1}^n y_i x_i - e^{a}\sum_{i=1}^n e^{x_i b} -\frac{b^2}{2\tau^2} \end{align*}Podemos considerar \begin{align*} h_1(a) = & a\sum_{i=1}^n y_i - e^{a}\sum_{i=1}^n e^{x_i b} -\frac{a^2}{2\sigma^2} \\ h_2(b)= & b\sum_{i=1}^n y_i x_i - e^{a}\sum_{i=1}^n e^{x_i b} -\frac{b^2}{2\tau^2} \end{align*} con lo que \begin{align*} h_1'(a) & = \sum_{i=1}^n y_i - e^{a}\sum_{i=1}^n e^{x_i b} -\frac{a}{\sigma^2} \\ h_2'(b) & = \sum_{i=1}^n y_i x_i - e^{a}\sum_{i=1}^n x_i e^{x_i b} -\frac{b}{\tau^2} \end{align*}
using RDatasets
prussian = dataset("pscl", "prussian")
┌ Info: Recompiling stale cache file /home/workstation/.julia/compiled/v1.2/RDatasets/JyIbx.ji for RDatasets [ce6b1742-4840-55fa-b093-852dadbb1d8b] └ @ Base loading.jl:1240
| Y | Year | Corp | |
|---|---|---|---|
| Int32 | Int32 | Cat… | |
| 1 | 0 | 75 | G |
| 2 | 2 | 76 | G |
| 3 | 2 | 77 | G |
| 4 | 1 | 78 | G |
| 5 | 0 | 79 | G |
| 6 | 0 | 80 | G |
| 7 | 1 | 81 | G |
| 8 | 1 | 82 | G |
| 9 | 0 | 83 | G |
| 10 | 3 | 84 | G |
| 11 | 0 | 85 | G |
| 12 | 2 | 86 | G |
| 13 | 1 | 87 | G |
| 14 | 0 | 88 | G |
| 15 | 0 | 89 | G |
| 16 | 1 | 90 | G |
| 17 | 0 | 91 | G |
| 18 | 1 | 92 | G |
| 19 | 0 | 93 | G |
| 20 | 1 | 94 | G |
| 21 | 0 | 75 | I |
| 22 | 0 | 76 | I |
| 23 | 0 | 77 | I |
| 24 | 2 | 78 | I |
| 25 | 0 | 79 | I |
| 26 | 3 | 80 | I |
| 27 | 0 | 81 | I |
| 28 | 2 | 82 | I |
| 29 | 0 | 83 | I |
| 30 | 0 | 84 | I |
| ⋮ | ⋮ | ⋮ | ⋮ |
# Año
x = sort(unique(prussian[!,:Year]))
20-element Array{Int32,1}:
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
# Número de muertes por patada de caballo
y = [ sum(prussian[!,:Y][ prussian[!,:Year] .== x[i] ]) for i in 1:length(x) ]
20-element Array{Int64,1}:
3
5
7
9
10
18
6
14
11
9
5
11
15
6
11
17
12
15
8
4
plot(p)
plot(p)
Sea $s\in\mathbb{R}$ y $X\sim \mathcal{L}\left( f(x) \right)$ con $f$ una función de densidad de probabilidad en $\mathbb{R}$ que es monótona en $[s,\infty)$.
Consideremos $U$ una variable aleatoria que toma valores en $(0,1]$, $\sigma >0$ y $Y=g(U)=s-\sigma\log(U)$ entonces $g^{-1}(y)=e^{ (s-y)/\sigma}$ y $U$ tiene densidad
\begin{align*} f_U(u) & = \frac{f_Y(y)}{\left| \frac{\mathrm{d}u}{\mathrm{d}y} \right|} = \frac{f_Y\left( y\right)}{\left| \frac{\mathrm{d}g^{-1}(y)}{\mathrm{d}y} \right|} = \frac{f_Y(y)}{\left| \frac{1}{g'\left( g^{-1}(y) \right) } \right|} = f_Y(y)\left| g'\left( g^{-1}(y) \right) \right| \\ & = f_Y(y) \left| -\frac{\sigma}{g^{-1}(y)} \right| =\sigma f_Y(y)e^{ (y-s)/ \sigma} = \frac{ \sigma f_Y(y) }{ u } \end{align*}Por construcción la variable aleatoria $Y$ toma valores en $[s,\infty )$. Por lo que podemos determinar la ley de probabilidad de $U$ al elegir $f_Y(y)=f$, nuestra densidad de probabilidad de interés. Así las cosas, si $(Z_1,Z_2)$ es tomado uniformemente debajo de la gráfica de $f_U(u)$ entonces $Z_1$ sería una realización de $\mathcal{L}(f)$ restringida a $[s,\infty )$, es decir usamos aceptación-rechazo. Este método es especialmente conveniente si con $f_Y=f$ se tiene que $f_U$ es estrictamente decreciente en $[s,\infty)$ ya que entonces podemos usar la cota $Mg(x)=f_U(s)$.
Si $f(x)=\frac{1}{\sqrt{2\pi}}e^{-x^2/2}$ entonces $f_U(u)=\frac{ \sigma }{\sqrt{2\pi}}e^{-y^2/2}e^{(y-s)/\sigma} $ es estrictamente decreciente si y sólo si
\begin{align*} \frac{\mathrm{d}}{\mathrm{d}y}\left( e^{-y^2/2 + (y-s)/\sigma} \right) <0 & \iff e^{-y^2/2 + (y-s)/\sigma}\left( -y + 1/\sigma \right) < 0 \\ & \iff \frac{1}{y} < \sigma. \end{align*}Recordemos que $s < y$ por lo que se puede elegir $\sigma\in [1/s, \infty )$. Así las cosas podemos utilizar $Mg(u)=\sigma e^{-s^2/2}/\sqrt{2\pi}$ y la regla de aceptación
\begin{align*} U_2 \leq \frac{f_U(U_1)}{f_U(s)}=\frac{ e^{-Y^2/2 + s^2/2} }{U_1}=e^{-Y^2/2 + s^2/2 + (Y-s)/\sigma}. \end{align*}La probabilidad de aceptación en este caso es un factor de $1/\sigma$ al ser $\sigma$ factor de $Mg$. Por esta razón tomamos el mínimo valor posible de $\sigma$ que es $1/s$.
1) Genera $U_1, U_2 \sim \text{Uniforme}(0,1)$ independientes y $Y=s - \frac{1}{s}\log(U_1)$.
2) Si $$U_2 \leq e^{-Y^2/2 + s^2/2 +s(Y-s) }$$ aceptamos $X=Y$, en caso contrario regresamos a 1).
Observando que $$U_2 \leq e^{-Y^2/2 + s^2/2 +s(Y-s) } \iff -2 log(U_2) \leq \left( \frac{\log(U_1)}{s} \right)^2$$ obtenemos un algoritmo menos costoso:
1) Genera $U_1, U_2 \sim \text{Uniforme}(0,1)$ independientes y $Y= -\frac{1}{s}\log(U_1)$.
2) Si $$-2\log(U_2) > Y^2$$ aceptamos $X=s+Y$, en caso contrario regresamos a 1).